function [corrplanes,sups,clusters]=a4C(data,e,u,r,q,minsup)
    [M,N]=size(data);
    corrplanes={};
    sups=[];
    Cores={};
    clusters={};
    for i=1:M
        o=data(i,:);
        o_neighbors=neighbors(data,o,e);
        [m,junk]=size(o_neighbors);
        if m<u
            continue;
        end
        [eigenvectors,eigenvalues]=PCA(o_neighbors);
        for t=1:length(eigenvalues)
            if eigenvalues(t) > q
                break;
            end
        end
        t=t-1;
        if t>=N-r
            MP=corrsimlarytymatrix(eigenvectors,eigenvalues,q);
            [count,junk]=size(Cores);
            Cores{count+1,1}=o;
            Cores{count+1,2}=MP;
            Cores{count+1,3}=o_neighbors;
        end
    end
    [CoreM,junk]=size(Cores);
    for i=1:CoreM
        O=Cores{i,1};
        O_MP=Cores{i,2};
        for j=i+1:CoreM
            P=Cores{j,1};
            P_MP=Cores{j,2};
            d1=sqrt((O-P)*O_MP*(O-P)');
            d2=sqrt((P-O)*P_MP*(P-O)');
            dist=max(d1,d2);
            if dist<e
                Cores{i,3}=[Cores{i,3};Cores{j,3}];
            end
        end
    end
    for i=1:CoreM
        points=Cores{i,3};
        [pm,junk]=size(points);
        if pm>=M*minsup
            [eigenvectors, eigenvalues]=PCA(points);
            cmean=mean(points);
            coe=eigenvectors(:,1)';
            const = -coe*cmean';
            plane=[coe,const];
            sup=pm/M;
            corrplanes{end+1}=plane;
            sups=[sups,sup];
            clusters{end+1}=points;
            unipoints=unique(points,'rows');
            unisup=length(unipoints)/M;
        end
    end
end

function [MP]=corrsimlarytymatrix(eigenvectors,eigenvalues,q)
    tk=50;
    for i=1:length(eigenvalues)
        if eigenvalues(i) > q
            eigenvalues(i)=1;
        else
            eigenvalues(i)=tk;
        end
    end
    MP=eigenvectors*diag(eigenvalues)*eigenvectors';
end

function [points]=MPneighbors(data,MP,core,e)
    
end

function [points]=neighbors(data, core, e)
    [M,N]=size(data);
    points=[];
    for i=1:M
        p=data(i,:);
        if isequal(p,core)
            continue;
        end
        if norm(core-p) <e
            points=[points;p];
        end
    end
end

% Find the smallest component vector
function [eigenvectors,eigenvalues]=PCA(data)
    [M,N]=size(data);
    mn=mean(data,1);
    newdata=data - repmat(mn,M,1);
    covariance=1/M*newdata'*newdata;
    [PC,V]=eig(covariance);
    V=diag(V);
    [eigenvalues, rindices]=sort(V);
    eigenvectors=zeros(N,N);
    for i=1:length(rindices)
        eigenvectors(:,i)=PC(:,rindices(i));
    end
end
